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ABSTRACT 

We use a model of the Galactic fountain to simulate the neutral-hydrogen emission 
of the Milky Way Galaxy. The model was developed to account for data on external 
galaxies with sensitive Hi data. For appropriate parameter values, the model repro- 
duces well the Hi emission observed at Intermediate Velocities. The optimal parame- 
ters imply that cool gas is ionised as it is blasted out of the disc, but becomes neutral 
when its vertical velocity has been reduced by 30 per cent. The parameters also 
imply that cooling of coronal gas in the wakes of fountain clouds transfers gas from 
the virial-temperature corona to the disc at ~ 2M0yr~^. This rate agrees, to within 
the uncertainties with the accretion rate required to sustain the Galaxy's star forma- 
tion without depleting the supply of interstellar gas. We predict the radial profile of 
accretion, which is an important input for models of Galactic chemical evolution. The 
parameter values required for the model to fit the Galaxy's Hi data are in excellent 
agreement with values estimated from external galaxies and hydrodynamical studies 
of cloud-corona interaction. Our model does not reproduce the observed H I emission 
at High Velocities, consistent with High Velocity Clouds being extragalactic in origin. 
If our model is correct, the structure of the Galaxy's outer Hi disc differs materially 
from that used previously to infer the distribution of dark matter on the Galaxy's 
outskirts. 
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1 INTRODUCTION 

Independent estimates of the baryon content of the Universe 
are provided by the theory of primordial nucleosynthesis 
(e.g. lPagell 19971 ) and studies of the inhomo Reneity of the cos - 
mic microwave background (CMB) fe.g. lSpergel et a"l]|2007^ . 
The resulting baryon density proves to exceed by a factor 
~ 10 that accounted for by the s tars and observed inter- 
stella r medium (ISM) of galaxies jProchaska &: TumlinsonI 
This dissonance suggests that most baryons are in 
intergalactic space. A fraction around 30% of these baryons 
may c omprise ionised gas associated with the local Lya 
forest (|Penton et al.l |2004| ). The rest is still missing, and 
probably comprises the warm-hot medium (WHIM) that 
is believed to permeate most of intergalactic space. In this 
picture, spiral galaxies, like more massive structures such 
as galaxy clusters, should be embedded in massive coro- 
nae of gas at the virial-temperature. These coronae should 
typically extend to a few hundred kiloparsecs from galaxy 
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centres l|Fukugita fc Peeblesll2006|V The obse rvational quest 
for this medium is ongoing (|Bregman|[2007l '). and there is 
some debate as to whether the observational constraints 
are compatible with the cosmological predictions or no t 
iRasmussen et al.ll2009l : lAnderson fc Bregmanlboid . I2OIII ). 

Cosmological coronae constitute a virtually infinite 
source of gas to feed star formation in galaxy discs. In 
fact several lines of evidence indicate that gas accretion 
from the intergalactic me dium (IGM) plays a n important 
role in galaxy evolution jSancisi et al.1 I2OO8I ). Studies of 
the stellar content of the Milky Way's disc show that 
the Star Formation Rate (SFR) in the solar neighborhood 
has declined b y a factor of onl y 2 — 3 ov e r the past 10 
Gyr jTwarogl [l980i: iRoch a-Pinto eTal] I2OOOI : ICignoni et al] 



l2008l : lAumer fc BinnevI [2009 1 . The slowness of the decline 



in the SFR suggests that the Galaxy's meagre stock of 
cold gas is constantly replenished. This suggestion is re- 
inforced by the scarcity of metal-poor G dwarfs, which 
aris es naturally if the Galaxy constantly accretes m etal-poor 
gas (|Pagel fc Patchettlll975hlChiappini et al.lll997l 'l. Finally, 
studies of the evolution of the cosmic SFR and thus the rate 
of gas consumption across the Hubble time show that at any 
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time galaxies must have accreted gas at a rate close to their 
SFR (|Hopkins etallbOOSl : iBauerineister et alJlioiol ) . 

While it is widely agreed that the WHIM contains 
the bulk of the baryons, and that disc galaxies sustain 
their star formation by accreting from the WHIM, there 
is no consensus at to how they accrete gas. Some au- 
thors have argued that cooling occurs spontaneously within 
WHIM coronae as a c onsequence of thermal instability 
l|Maller fc Bullock |2004| ). producin g cold clouds similar to 
the G alactic high velocity clouds (|Wakker fc van WoerdenI 
Il99'i1 ). However, linear perturbation analysis shows that in 
a corona stratified in a galactic potential, buoyancy can 
very efficiently suppress th ese instabil ities jMalagoli et alj 
ll987l : lBalbus fc Sokeij[l989l ). Recentlv. [Binnev et all (|2009l ) 
employed Malagoli et al.'s method to study the growth of 
thermal instabilities in cosmological coronae typical of spiral 
galaxies with different masses. They found that in each case, 
the combination of buoyancy and thermal conduct ion to- 
tally s uppresses the growth of thermal perturbations. iNipotil 
l|201Cl ) studied the growth of such perturbations in a more 
realistic rotating corona, and re ached similar concl usions. 
These results were confirmed by Ijoung et all (|201ll ) using 
hydrodynamical simulations based on adaptive mesh refine- 
ment. They are also broadly consistent with the findings 
of the numerical simulations presented bv lKaufmann et al.l 
which show that a corona must have a peculiarly 
flat profile of specific entropy to produce cold clouds. How- 
ever, Kaufmann et a l.'s simulations did not include ther- 
mal conduction, and iBinnev et all (|2009l ') found that con- 
duction suppresses perturbations in flat-entropy coronae. It 
appears therefore rather unlikely that instabilities can grow 
anywhere in WHIM coronae to form clouds that will even- 
tually feed the star formation in the disc. 

If coronae were thermally unstable, one would expect 
to find them studded with clouds of cool gas. Searches 
of ever-increasing sensitivity have been undertaken for HI 
clouds around external galaxies. These searches have con- 
sistently failed to discover massive HI cloud s at large dis- 
tance s from galaxies, eve n in gal axy groups (ILo fc Sargent! 
1 19791 : [Pisano ot al. 20071 : llrwin et aL 20091'). Deep surveys 
in HI emission such as HIPASS and ALFALFA, indicate 
that massive bodies of H I are always ass ociated with stars 
and thus constitute t he ISM of a galaxy (|Dovle et al.ll2005l : 
ISaintonge et al.ll2008l ). 

What sensitive observations of nearby galaxies have 
revealed is that star-forming galaxies hold significant 
quant ities of HI one or more k i loparsecs above their 
discs (Ivan der Hulst fc Sancisil [l98l : iFraternah et al.ll2002l : 
lOosterloo et al.l2007 ). This gas rotates around the centre of 
the galaxy but at lower speed than the gas in the disc, and 
there are indications that it ha s a net component of velocity 
inwards (|Fraternah et al.ll2002l : iBarbieri et aLlbOOSl l. 

Since the temperature of HI has to be lower than the 
virial temperature by at least a factor 100, this HI cannot 
be kept off the disc by pressure forces, but must be mov- 
ing essentially ballistically. Numerical simulations of the im- 
pact of supernovae on galactic discs predict that most of the 
gas ejected from the disc by a superbubble s hould indeed 
be m uch cooler than the virial temperature (|Melioli et al.l 
l2008h . Observationally, there are strong indications that in 
both the Milky Way and external galaxies HI has been 
blasted out of the disc by supernova-powered superbubbles 



- prominent extraplanar HI clouds occur near regions of 
the disc in which the surface densi ty of HI is depressed 
and there are m assive young stars l|Boomsma et al.| l2008l : 
IPidoprvhora et al...2007i ). Another indication that star for- 
mation gives rise to extraplanar HI is the observation that 
galaxies such as NGC 891 with unusually large SFRs also 
hold anomalously large fractions of their H I as extraplanar 
gas. Consequently, there is strong evidence that extraplanar 
H I is gas that has been blasted out of the disc by supernova- 
powered superbubbles, and is orbiting in the galaxy's gravi- 
tational field. That is, extraplanar H I constitutes a galactic 
fountain. 

iFraternali fc BinnevI (|2006l . |2008| . hereafter FB06, 
FB08) modelled galactic fountains with and without accre- 
tion from the surrounding environment, and applied their 
models to NGC 891 and NGC 2403. They found that the 
kinematics of the extra-planar gas - rotational lag with re- 
spect to the disc and global inflow - could only be explained 
once accretion from the corona was included. Moreover, the 
accretion rate required to explain the observed kinematics 
turned out to be very close to the SFR in these galaxies. 
An aspect of this work that did not make sense physically 
was the assumption that as a cloud of HI moved through 
the virial-temperature corona, its mass had to increase as a 
result of accreting coronal gas. Since the cooling time of the 
corona greatly exceeds the time required for each parcel of 
coronal gas to flow past the cloud, accretion should not take 
place^ 

iMarinacci et al.l (|2010l ) investigated the interaction be- 
tween the H I clouds of the fountain and the ambient corona 
using 2D hydrodynamical simulations. They found that as 
the cloud moves, gas is stripped off it and mixed with the 
much hotter coronal gas. If the metallicity and pressure of 
the mixed gas is high enough, radiative cooling is fast and 
the mixture cools to ~ 10* K within a dynamical time. Then 
the gradually eroding H I cloud trails a wake of cold gas be- 
hind it, and the total mass of HI increases with time , just 
as the models of FB08 required. IMarinacci et al.l (|2010l ) sug- 
gested that this is the mechanism by which gas is transferred 
from the WHIM to star-forming discs. 

Our location within the disc of the Milky Way affords us 
a radically different view of the local disc/corona interaction 
from views we have of such interactions in external galaxies. 
Extraplanar gas was first observed around the Milky Way i n 
the form of High Velocity Clouds (HVCs) (e.g. IOortllT966l l. 
For decades after the discovery of HVCs, t heir distances and 
hence masses remained controversial (e.g. iBlitz et al.lll999l : 
Ide Heii et al.l [2002). In the last decade it has been estab- 
lished that the original HVCs are either associated with 
other galaxies (the Magellanic Clouds or M31) or they lie 
within about 10 kp c from the Galactic disc an d consequently 
have small masses (|Wakker et al.ll2007l . [2003 ). However, the 
HVCs represent only the tip of the iceberg of extraplanar 
gas: most of this gas should be detected at Intermediate 
Velocities, an d indeed at such velocities there is abundant 
HI emission (jMarasco fc Fraternaiill201ll . hereafter MFll). 
Part of this emission constitutes the Intermediate Velocity 
Clouds (IVCs), bright complexes with disc-like metallicity 
located around ~ 2 kpc from the Sun, whi ch have been r e- 
garded as nearby Galactic fountain clouds (|Wakkejl200lh . 

MFll defined a kinematic model of the Galaxy's HI 
layer and used it to simulate the HI datacube of the 
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Leid en-Argentine-Bonn (L AB) survey of Galactic H I emis- 
sion l|Kalberla et al.ll2005l ). This model reproduced well the 
Galactic emission at intermediate velocities, including the 
IVCs, which are interpreted as a local manifestation of a 
more global phenomenon. In this paper we ask whether we 
can obtain better-fitting datacubes by replacing the kine- 
matic model of MFll with the dynamical model that FB06 
and FB08 developed to account for observations of external 
galaxies. In Section[2]we describe the model we will be using, 
which is an updated version of that described in FB08, and 
explain how we evaluate the fit between a model datacube 
and that of the LAB survey. In Section 13.11 we investigate 
whether the H I thick disc of the Milky Way can be produced 
by a pure galactic fountain. Then in Section [3.21 we present 
evidence that this fountain interacts with the ambient coro- 
nal gas. In Section|4]we use our model to infer the large-scale 
structure of the Galaxy's H I halo, and argue that there is a 
fundamental distinction between HVCs, which have an ex- 
tragalactic origin, and Intermediate Velocity Clouds, which 
are dominated by fountain gas. In Section[5]we show that our 
model's parameter values are remarkably consistent with re- 
sults from earlier studies of (a) small-scale hydrodynamical 
simulations of turbulent mixing, and (b) analysis of data for 
external galaxies. Section[5]sums up and looks to the future. 



2 THE MODEL 

We set up a galactic fountain model for the Milky Way 
by integrating orbits of gas clouds in the Galactic poten- 
tial. The latter has been constructed using the standard de- 
composition into a bulge, stellar and gaseous discs, and a 
dark-matter halo. The physical parameters and functional 
forms for the various com ponents have been taken from 
iBinnev fc Tremaind l|2008l . Model II, p. 113). 

The clouds are ejected from the disc into the halo, and 
we follow their orbits until they return to the disc as de- 
scribed in FB06. The probability of ejection at speed v at 
angle 6 with respect to the vertical direction is 

^(^'^^"'^^^(-s^^fi^)' « 

where is a characteristic velocity and F is a constant that 
determines the extent to which clouds are ejected perpen- 
dicular to the disc. Here we treat hv as a free parameter but 
adopt r = 5 from FB06 - such a large value of F implies 
strong coUimation of ejecta towards the normal to the plane 
by cool gas near the plane. The integration is performed in 
the {R, z) plane and then spread through the Galaxy by as- 
suming azimuthal symmetry. The number of clouds ejected 
as a function of radius is proportional to the SFR density at 
that radius (see Section [2]2}. Other parameters and proper- 
ties of the model are as described in FB06 except for some 
updates to the model that are specified in the following sub- 
sections. 

2.1 Phase-change 

We follow the orbits of gas that leaves the disc at relatively 
low temperatures ( ^ 10"* K) after being swept up in the ex- 
pansion of a superbubble that is powered by multiple su- 
pernova explosions. The virial-temperature gas filling the 



bubble is presumed to contain negligible mass and to merge 
with the pre-existing coronal gas. In the Milky Way, emis- 
sion from the H I halo occurs mainly at negative line-of-sight 
velocities (|van Woerden et al.lll985l ). This observation sug- 
gests that the neutral gas is seen mainly as it descends to 
the plane, so a significant fraction of the ejected gas must 
be ionised, and indeed from their 3D kinematic models of 
the galactic halo MFll estimated that the rising clouds are 
~ 60 per cent ionised. If we are to compare our models with 
the LAB survey of neutral gas, we have to find a way to 
estimate the fraction of a cloud that is ionised at each part 
of its orbit. 

FB06 built models with and without phase-change, the 
latter having the whole orbit "visible" as H I gas, the former 
only the descending part. Here we refine this treatment as 
follows. We assume that a cloud ejected from the disc with 
a kick velocity Ukick will become visible (i.e. neutral) only 
when 

^^z(t)<W.,0(l-/ion), (2) 

where is the vertical component of the cloud's velocity, 
Vz,o ~ Wkick X COS 6 and ^ /ion ^ 1 is a free parameter 
that regulates the "visibility" of the cloud. If /ion = 1, the 
fountain cloud is neutral only for negative values of Vz (i.e. 
in the descending part of the orbit), while if /ion = 0, the 
cloud will always be visible. We infer the correct value of 
/ion from the data cube (see Sections |3 . 1 1 and . 

2.2 Star- Formation Law 

We assume that the strength of the supernova feedback, i.e. 
the number of gas clouds ejected per surface element, is pro- 
portional to the local SFR. Given that our model is axisym- 
metric, we require an estimate of the SFR as a function of 
R. In FB06 and FB08 that was estimated from the surface 
density of neut ral plus molecular gas us i ng th e Schmidt- 
Kennicutt law (|Schmidtl Il959l : lKennicutl ll998^ . The SFR 
was then set to zero beyond a certain cut-off radius. Here, 
we refine this procedure as follows. We use a sample of galax- 
ies with known radial trends of gas surface density and SFR 
density and derive a star-formation law for molecular gas. 
We then use this law to obtain the Milky Way's SFR as a 
function of R from th e axisymmetrised su r face d ensity of 
molecular gas given in iBinnev fc Merrifieldl l|l998l ). 

Fig. [1] is a plot of SFR density versus the molecular gas 
density. The fi lled triangles show data points for galaxies in 
the sample of iLerov et al.l (|2008l ) that have stellar masses 
above 10^ Mq. Each point refers to an azimuthal average at 
the given radius. The scatter in the relation between gas 
density and SFR is remarkably small considering that the 
points come from very different galaxies. We fitted these 
points with a power law of the form: 

where Espr and are, respectively, the SFR and the 

molecular gas density. We find 

A = 10.8 ± 2.3 X 10"" Mq yr"^ kpc"^ 

TV = 0.85 ±0.04. (4) 

The slope of this relation is very different from the stan- 
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Figure 1. Star formation rate densities vs gas surface densi- 
ties, both for molecular gas (filled triangles) and neutral atomic 
gas (crosses); fro m the sample of nearby galaxies studied by 
iLerov et al. 1 1|2008| '). The points come from an azimuthal averages 
at different radii over the galactic discs of 17 galaxies. The curve 
shows the power law fit described in the text. 



dard Schmidt-Kennicutt law, mainl y because we are con - 
sidering only molecular gas (see also Krumholz et al.ll201ll ). 
Including HI would increase the slope to about 1.3 — 1.4, 
however the average H I surface density correlates very little 
with the SFR d e nsity - see the grey crosses in Fig. [1] and 
iKennicufrt et al.1 l|2007l ). Therefore, given that the distribu- 
tion of molecular gas in the disc of the Milky Way is known 
reliably, we prefer to obtain the SFR density from equation 
([3]). Our results do not depend strongly on the SF law. 



2.3 Supernova-driven gas accretion 

The pure galactic fountain model described above is modi- 
fied to include interaction with the coronal gas according to 
FB08. We as sume that the m a in me chanism at work is that 
described by iMarinacci et al] l|2010l ): the Kelvin-Helmholtz 
instability generates a turbulent wake behind each cloud in 
which material stripped from the cloud mixes with coronal 
gas, enhancing the metallicity and density of the latter and 
substantially decreasing its cooling time. As a consequence, 
part of the corona condenses into the wake and the mass of 
"cold" gas increases. The combined mass of cold gas in the 
cloud and the wake grows with distance along the cloud's 
path through the halo very much as proposed by FB08. We 
demonstrate this quantitatively in Section [5. II 

Condensation in the wake modifies the kinematics of 
the cold gas in a way that depends on the kinematics of 
the corona. The latter is not known a priori, but must be 
strongly infiuenced by this interaction because the portion of 
the corona that interacts with fountain clouds does not con- 
tain much mass, and the interaction has been taking place 
throughout the disc's lifetime. FB08 considered the effect of 
the drag between the clouds and the corona and concluded 
that, if there were no exchange of gas between the clouds 
and the corona, ram pressure between the cloud and the 
corona would force the corona to corotation with the cold 
gas in a time shorter than the fountain's dynamical time. 



In reality much of the coronal gas that acquires momentum 
from a fountain cloud subsequently condenses in the cloud's 
wake, so the momentum that is transferred to the coronal 
gas is in the end not lost by the whole body of cool gas. 
IMarinacci et al.l (|201ll ') studied the transfer of momentum 
between fountain clouds and the corona in hydrodynamical 
simulations and found that there is a net transfer of momen- 
tum to the corona until the speed at which the clouds move 
through the corona falls to 50 — lOOkms"^, depending on 
the physical properties of the system. Below this threshold 
speed, net momentum transfer to the corona ceases as a re- 
sult of coronal gas cooling in the cloud's wake. In view of 
the short timescale on which ram-pressure can change the 
rotation rate of the corona, it is natural to assume that the 
corona has reached the velocity at which all the angular mo- 
mentum imparted by ram pressure is subsequently recovered 
through condensation of coronal gas in the wake. Thus we 
assume that the rotation velocity of the corona is 



v^{R, z) — \ R 



OR 



(5) 



where ^{R, z) is the potential of the Galaxy and iiiag = 
75kms~^ is the (constant) velocity offset between clouds 
and corona at which condensation recaptures angular mo- 
mentum transferred by ram pressure. 

In this picture, the mass of cold gas (cloud-|-wake) as- 
sociated with any cloud increases as 



am. 



(6) 



where ct is the condensation rate, while the combined effect 
of condensation and drag decreases the velocity centroid of 
the cold gas at a rate (FB08) 



^drag (1 H" ^/^drag) 



{v - v^) 



(7) 



where tdrag is the time at which the relative velocity between 
the cold gas and the corona halves due to the drag only. Since 



"lag' 



the drag efficiency decreases as the 



kinematics of the corona approaches that of the disc. We 
assumed a c onstant value for tdrag o f 800 Myr, according to 
the results of lMarinacci et ahl (|201ll ) for |u — ii^l = 75kms~^ 
(see Section r5.1|) . and we fitted the parameter a. As we see 
below (Section 13. 2p . in order to reproduce the data we need 
a > 1/tdrag, i.e. the condensation is the dominant process 
at work, while the drag plays only a secondary role. Note 
that fixing tdrag at 800 Myr implies that we are considering 
clouds with masses and sizes similar to those of Marinacci 
et al.'s simulations (~ 10* Mq, 100 pc). 



2.4 Comparison with the data 

The comparison between the models and the data is per- 
formed by building model datacubes that resemble the LAB 
datacube (see FB06 and MFll). After each timestep of the 
integration, the positions and velocities of the particles are 
projected along the line of sight of an observer located in 
the mid-plane of the Galaxy, at a radius Rq = 8.5 kpc, and 
moving at a circular speed ulsr = 220 km s^'^. Using the 
same gas density for each cloud, we construct from these 
positions and velocities an artificial datacube. The artificial 
datacube is smoothed to 8° resolution and then compared 
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Figure 2. Longitude- velocity diagrams at latitudes ±30° above and below the mid-plane of the Galaxy. Left panels: LAB survey. Middle 
panel: disc model used throughout this paper. Right panel: LAB survey after the region of the disc has been removed, the remaining 
emission is used to compared our models with the data. Each datacube is smoothed to 8° resolution. Contour levels in brightness 
temperature range from 0.04K to 81.92 K scaling by a factor 2. 



with the LAB datacube at the same resolution. The compar- 
ison is performed only in the portion of (I, b, wlos) space to 
which the thin disc should not contribute. Following I Wakken 
l|l99ll ), we use the deviation velocity wdev to identify this 
region. We model the footprint of the thin disc as described 
in MFll and exclude data with |fDEv| < 35kms~^. We 
also excluded the H I emission of the Magellanic Clouds and 
Stream, the Stream's Leading Arm, the GCP complex, the 
Outer Arm and all external galaxies. In the following we 
refer to the surviving part of hlos) space as the halo 
region. The total model flux in the halo region is normalized 
to the corresponding LAB flux. This step fixes the mass of 
the HI halo. 

Although the quantitative comparison between mod- 
els and data is done only in the halo region, in all the 
plots of this paper a thin disc has been added for presen- 
tation purposes. Thi s disc has a density profile taken from 
Binney fc Merrifieldl (|l998l ) and a scale-height taken from 
Kalbcrla ct al. (20071). The velocity dispersion decreases lin- 
early from 12kms~^ in the Galactic centre to 6kms~^ at 
Ro, and it remains constant at larger radii. 

The leftmost column of Fig. [2] shows two longitude- 
velocity {l,v) plots at latitudes ±30° extracted from the 
LAB survey. The middle column of the same figure shows 
our model disc and the rightmost column shows the emis- 
sion in the halo region of the LAB datacube. All the plots 
are centered at the anti-centre (i = 180°). Note that our cut 
is rather conservative and there should be no contamination 
from disc emission in the halo region. 

In the application of the model to external galaxies, the 
fountain clouds have been modelled as point particles be- 
cause individual clouds, with typical sizes ~ 100 pc, were 
unresolved in the data (FB06). In the Milky Way nearby 



clouds and their wakes will be resolved, so we have to check 
at each time step whether the combined emission of a cloud 
and its wake is resolved at the angular resolution of the 
data. If the emission is resolved, we spread the flux in a cir- 
cular region around its centre. The flux density is assumed 
to decrease Gaussianly with distance from the centre, with 
FWHM equal to the apparent size of the cloud and its wake. 
Given that the data are used at 8° resolution, a cloud and 
its wake are resolved only if their distance is smaller than 
700 X pc, where Dd i s the system's lengthsc ale. Ac- 

cording to the simulations of lMarinacci et al.1 l|2010t ). turbu- 
lent wakes have sizes of 1 — 2kpc, so they will be resolved 
out to distances of about lOkpc. The pattern of emission 
from a cloud and its wake is quite elongated, having axis 
ratio '-^ 0.5, but the long axis of the emission is randomly 
orientated, so our circular Gaussian smoothing should pro- 
vide an adequate representation of the total emission from 
an ensemble of clouds. We tested the effect of resolution by 
increasing Del up to 2kpc and found that the model dat- 
acube looked just like a smoothed version of that produced 
with smaller values of Dd- In the following we use models 
with small clouds, Del ~ 200 pc, for which the datacube can 
be computed signiflcantly more quickly. 

Once we have a model datacubes at the same res- 
olutions and with the same total flux as the LAB dat- 
acube, we compare them quantitatively by calculating the 
residuals between them. This calculation is performed by 
adding up the differences between the models and the data, 
pixel by pixel, in the halo region. We used absolute differ- 
ences, squared differences and also weighted residuals, i.e. 
data— model|/(data±model). The residuals are calculated 
for a set of models with different input parameters such as 
the kick velocities of the fountain clouds, and the model 



6 A. Marasco, F. Fraternali, & J. J. Binney 



Idata-modI 
Idata-modl^ 



ldata-modl/{data+mod) 



to 

■g 

to 

w 

"O 

tu 
to 

16 
E 
o 



1.16 
1.14 
1.12 
1.1 
1.08 
1.06 
1.04 - 
1.02 - 
1 - 
0.98 - 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
'ion 

Figure 3. Resitiuals between motiels anil data as a function of 
/ion for the pure galactic fountain with hv = 70kms~^, eval- 
uated using absolute differences (solid line), squared differences 
(short-dashed line) and weighted differences (long-dashed lines). 
All residuals have been divided by their respective minimum. 



with the smallest residuals is our "best model". In the fol- 
lowing section we describe our best pure fountain model and 
our best fountain + condensation model, showing that the 
latter provides a better description of the Galaxy's H I data. 



3 RESULTS 

3.1 A pure galactic fountain 

In a pure galactic fountain, the angular momentum of clouds 
is constant and orbits depend only on the characteristic kick 
velocity hv of the clouds. However, the shape of the model 
datacube is also affected by the ionised fraction /ion, which 
has to be treated as a freely variable parameter. 

We have calculated the residuals for in the range 
(30, 120) kms~^ and /ion between and 1. The weighted and 
non-weighted residuals give different results but the confi- 
dence contours overlap for hv — 70kms~^ and /ion = 1- 
In Fig. [3] we show the normalized residuals as a function 
of /ion for /iv — 70kms~^. Clearly the differences between 
models and data decrease as /ion approaches 1, regardless of 
the method of comparison. The mass of the H 1 halo derived 
for this model is A/haio ~ 5.2 x 1O*A^0. Here and below the 
halo masses are estimated by excluding fountain clouds at 
\z\ < 0.4kpc, which is roughly th e blowout height of a su- 
perbubble (see lSpitoni et al]|2008h . Taking into account also 
clouds at heights 0.2< \z\ <0.4kpc, the Hi halo mass would 
become Mhaio = 7.2 x 10* M©. 

Fig. |4] shows two representative l,v diagrams (centred 
on the anticentre) at latitudes b = 60° and 6 = 30° for: the 
best pure galactic fountain model (middle column); the same 
model without phase-change (left column); the data (right 
column). The difference between the two models is striking, 
especially at higher latitudes. The l,v plots for the neutral 
(./ion ~ 0) fountain are rather symmetric with respect to 
the zero velocity line. In the ionised fountain by contrast, 
the gas appears systematically located at negative veloci- 
ties. This effect is a consequence of clouds of the ionised 



fountain being visible only as they fall back to the disc. The 
data clearly display this preference for more negative than 
positive velocities. 

The leftmost column in Fig. [6] shows additional l,v di- 
agrams for the best pure fountain model. It reproduces the 
most prominent features, so it gives a good description of the 
data. It is interesting to compare Fig.[6]with Fig. 6 of MFll, 
which shows the prediction of an optimised kinematic model. 
In several locations the dynamical model fits the data better 
than the best kinematic model of MFllQ For example, at 
latitudes b — ±15° the superiority of the present model is 
very clear in the central regions of the Galaxy, 0° < I < 60° 
and 300° < I < 360°. At latitudes b = ±30°, the superiority 
is clear at / ~ 270° (-90° in MFll). 

The kinematic model was built under the assumption 
that the scaleheight and kinematic parameters of the halo 
(rotation velocity gradient, vertical and radial velocities) are 
independent of radius. Although the assumption is restric- 
tive, the model has been fitted to the data by freely varying 
these quantities, regardless of their dynamical plausibility. 
The fact that a dynamical model with only two free parame- 
ters fits the data better is remarkable, and strongly supports 
the idea that the H 1 halo of the Milky Way is produced by 
supernova feedback from the Galactic disc. 



3.2 Including the condensation of the corona 

The best pure fountain model described above is extreme 
in that it requires the full ionisation of ejected clouds and 
a large mass of the gaseous halo (see Table [1]). Moreover 
there are several locations in the datacube where the model 
systematically fails to reproduce the data in detail. The ar- 
rows in Fig. [6] indicate such regions. In these regions radial 
motions are important. For instance filled arrows show the 
region around the anticentre, 150° < I < 250°, where the 
data, almost at every latitude, have a broader distribution 
in velocity than the pure fountain model. Gas in this region 
is flowing towards the centre of the Galaxy from beyond the 
Solar circle. In a pure galactic fountain, orbits seldom move 
in this way (FB06) and independently of the ionised fraction 
they do not produce much emission in the marked region. 

We now show that including the interaction be- 
tween clouds and the corona as described in Section 12.31 
substantially improves this situation. The interaction is 
parametrised by the specific condensation rate a. The larger 
a is, the faster gas condenses from the corona, the greater 
the rotational lag of the Hi, and the more prominent is in- 
ward radial motion of Hi (see FB08). 

Now the residuals between the data and the models are 
minimised by varying three parameters (/iv, /ion, a). We as- 
sume that condensation and drag are suppressed close to the 
Galactic plane (|2| < 0.4 kpc). A full analysis of this three- 
dimensional parameter space would be too expensive com- 
putationally, so we explore only the space (a, /ion) and fix the 
kick velocities at three values: = 60, 70 and 80kms~^. 
For each value of we build a map of the residuals be- 
tween models and data by varying the other parameters in 
the ranges < /ion < 1 and -3 < log(a/Gyr"^) < -1.67 
with a step of 0.033 in both directions. This map is then 



^ Note that the plots in MFll are centred at the galactic centre. 
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Figure 4. Longitude- velocity (l,v) diagrams at 6 = 60° (top panels) and fe = 30° (bottom panels). First column: pure fountain model 
without phase-change (/jon = 0); second column: best pure fountain model (/ion = 1); third column: the LAB data. Each datacubc is 
smoothed to 8° resolution. Contour levels in brightness temperature range from 0.04 K to 81.92 K scaling by a factor 2. 



Table 1. Best galactic fountain models for the Milky Way. 



Fountain model 






^^^HI+HII 


/ion 


a 




Residuals'^ 




(kms-l) 








(Gyr-i) 


(Moyr-i) 




Neutral (pure) 


70'' 


4.6 X 10* 


4.6 X 10* 


0.0* 


O.O'' 


0.0 


1.60 


Ionised (pure) 


70 


5.2 X 10* 


10.4 X 10* 


1.0 


0.0'' 


0.0 


1.48 


With corona condensation 


70 


2.7 X 10* 


3.0 X 10* 


0.3 


6.3 


1.6 


1.26 



" estimated above 400 pc from the plane; '' fixed; calculated as absolute differences, the closer to 1 the 
better the fit. 
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Figure 5. Contour plots of residuals between models with corona condensation and data evaluated in the parameter space (log(Q),/ion) 
for hv = 70kms~^. Different panels compare models and data in different ways, as labelled on bottom. In each panel all values have 
been divided by the respective minimum. The fourth panel is the sum of the previous three, divided by the resulting minimum. Contour 
levels at 1.005, 1.01, 1.02, 1.05. 
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Figure 6. I, v slices at 6 different latitudes indicated at the top right corner of the rightmost plots. First column: best pure Galactic 
fountain; second column: best fountain with condensation of the corona; third column: the LAB data. Each datacube is smoothed to 8° 
resolution. Contour levels in brightness temperature range from 0.04 K to 81.92 K scaling by a factor 2. The arrows show regions where 
radial motions are present and the model with coronal condensation reproduces the data better than the pure fountain model, see text. 
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smoothed to the resolution of 0.1 x 0.1 in order to min- 
imise the stochastic effects due to the probability distribu- 
tion function (eq.[T]) used to build our models. As in the case 
of a pure fountain, different residuals give different results. 
Regardless of a and /ion, absolute and squared differences 
have lower values for — 80kms~^, while relative residu- 
als favour hv = 60kms~^. Hence we set our kick velocity in 
the middle and focus on models with hv = 70kms~^. 

The first three panels in figure[5]show the residual maps 
obtained with — 70kms~^ for different types of residu- 
als. In each panel all values are divided by the respective 
minimum, thus they are dimensionless and ^ 1. Combining 
these three panels we obtain the fourth (rightmost) panel, 
which is the sum of the previous three divided by the result- 
ing minimum, which occurs at /ion = 0.3 and log(Q:) — —2.2 
(so Q = 6.3Gyr~^). For this best model the mass of the HI 
halo is 2.7 x 10* M©. Table [1] lists the parameters. 

The middle column of Fig. [6] shows that the model with 
coronal condensation reproduces the emission near the anti- 
centre better than the pure fountain model. The open arrows 
(latitudes b = ±30°) indicate other regions where the inclu- 
sion of coronal condensation improves the fit to the data. 
These are regions in which H I is flowing out from the inner 
disc, and the pure fountain provides too much emission. 

For our best-fitting value of a, the condensation rate 
of coronal gas into the clouds' wakes is M = 1.6M0yr~^. 
Correcting this value for the He content increases the con- 
densation rate to 2.3 Mq yT~^ . This value is remarkably close 
to_the_Galaxy's SFR, which we assume to be ~ 3Mq yr~^ 
jPiehl et all 12006). of which ~ IMgyr"^ should be ac- 
counted for by gas ejected from stars in the disc. We con- 
sider it highly significant that M deduced from an H I survey 
should agree so well with that required to sustain the SFR, 
which is obtained without reference to HI. 

The exact value of the condensation rate M depends 
on several factors. An important parameter is the value of 
flag (eq. El, which regulates the rotation of the corona and 
therefore the effici ency of both the drag (tdrag oc v^^) and 
the condensation (Marinac ci et al.ll201ll ). However, experi- 
ments show that the dependence of M on the adopted value 
of Wiag is not strong: (uiag — 50kms~^, tdrag = 1200 Myr) 
yields M = 2.8M0yr"^ (3.9M0yr"^ when corrected for 
He), while (i^iag ~ lOOkms^^, tdrag ~ 600Myr) yields 
M = 1.0 Mq yr"^ (1.4 Mq yr"^ when corrected for He). It is 
interesting to note that M depends very little on the crite- 
rion used to minimise the residuals. Although the three cri- 
teria yield rather different values for the best-fitting a (see 
Fig. [5|, M only varies between 1.5 and 1.7M0yr~^ given 
that the halo masses also vary between the three models. 

Our parametrisation assumes that a is spatially con- 
stant. Fortunately, the density of the corona should not vary 
too much in the regi ons where th e interaction with the foun- 
tain is efficient t Marinacci et al. l 2011), so spatial variation 
of a is not expected to give significantly different results. 

On the whole, the dynamical modelling of the Galactic 
HI halo performed in this work is a success. The results 
are consistent with the conclusions MFll drew from their 
kinematic models. They found negative values for both the 
average vertical and the radial velocities of the halo material, 
and deduced that ejected clouds must be partially ionised 
(but in a fraction larger than our current determination). 
Our dynamical analysis identifies these properties as arising 
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Figure 7. Rotational velocities for the H I layer of the Milky Way 
at different heights from the plane (thin lines), compared to the 
disc rotation curve (thick line). The velocities are derived from 
our best mixing fountain model by taking the weighted average 
of at that {R, z) location. 



from the sharing of angular momentum with gas accreted 
from the lagging corona. Finally, our Galactic fountain is 
dynamically sustainable. It requires a rate of kinetic-energy 
injection equal to 6.2 x lO'^^ergs"^, which with a supernova 
rate SNR= O.O3M0yr~^ and an energy per supernova of 
1 X lO^^erg corresponds to an efficiency of about 0.7%. 



4 PROPERTIES OF THE GALACTIC H I 
LAYER 

Now that we have a model that reproduces all the main 
features of the LAB datacube, we can investigate in detail 
the physical properties of the Galaxy's extraplanar H I layer. 

4.1 Rotation versus height 

Fig.[7]shows the rotation curve of the Galaxy's extra-planar 
gas at different distances from the plane as a function of 
R. These values were derived from our best model by sim- 
ply taking the weighted mean of the azimuthal components 
of the velocities of the particles at different heights. As ex- 
pected, the halo gas is lagging the rotation of the disc (black 
curve). Roughly 50% of this lag is due to gravity (FB06) 
whilst the other half is produced by a combination of the 
interaction with the corona and the phase-change. At Ro 
the gradient is about —14.3 ± 1.1 kms"'^ kpc~^ in excellent 
agreement with the average value derived in the inner re- 
gions of the Milky Way by MFll (-15 ± 4kms"^ kpc"^), 
which also coinc i des w ith the value derived for NGC 891 by 
lOosterloo eTsI] l|2007l ). 

4.2 Thickness of the HI layer 

A second fundamental property of the Galaxy's HI halo is 
its thickness. To obtain it we fitted the vertical density pro- 
files in our models at different radii with exponential func- 
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Figure 8. Scalo-hoight as a function of R for our best halo model 
including condensation of the coronal gas. 



tions. Fig. [8] shows the (exponential) scale-height of the HI 
halo of our best model as a function of R. The thickness 
increases with R because the gravitational restoring force 
to the plane diminishes outwards (FB06). The halo density 
decreases abruptly for 7? > 14 kpc, so its thickness cannot be 
reliably determined at larger radii. The shape of the scale- 
height as a function of R is partially due to the assumption 
of our model that the kick velocity does not change with 
radius (for a discussion see FB06). In the inner Galaxy there 
has long been evidence for a populatio n of H I clouds extend- 
ing up to ^ 1 kpc above the midplane (iLockman 1984, 2002; 



iFord et al.) I2OI0I '). iMarasco fc Fraternaliri|201ll ) assumed a 
constant thickness for the halo and derived a value of 1.6 kpc 
using a sech'^ formula, which corresponds to ~ 800 pc for an 
exponential function. This value agrees well with the average 
of t he scaleheights p lotted in Fig. [8] 

iKalberla et all l|2007l ) studied the vertical structure of 
the Galaxy's HI layer on the assumption that the layer is 
in hydrostatic equilibrium. Using this assumption, they in- 
ferred that while at Rq the HI layer extends only a few 
hundred parsecs from the plane, at i? ~ 35 kpc the verti- 
cal structure of the HI layer can be fitted by a Gaussian 
distribution with scale height 2.5 kpc. Thus they deduced 
a very extended H I layer with strong flaring. They also in- 
ferred the existence of a massive ring of dark matter between 
7? = 13 and 18.5 kpc. The assumption of hydrostatic equilib- 
rium leads, h owever, to models that fail to match the d ata 
correctly (see Barnabe et ai]|2006l : iMarinacci et al.ll2010l ). If 
[Kal bcrla et al. (2007) had used the prediction of a fountain 
model for the distribution of H I near _Ro rather than a model 
based on hydrostatic equilibrium, they would have come to 
different conclusions regarding the radii responsible for each 
quantity of emission in the HI datacube. Consequently, the 
distribution of the matter and the structure of the H I flare in 
the Milky Way should be re-derived. While this topic merits 
further work, our comparisons of the fountain model with 
the data yield no compelling evidence for substantial flaring 
of the HI layer much beyond R ^ 15 kpc. 





1x10"^ - 


C\J 


8x10^ 


kpc" 


6x10"'^ 




4x10-^ - 







2x10"* - 




0x10""° - 




-2x10"^ 




1x10'^ - 




8x10'^ - 


kpc" 


6x10"^ - 




4x10'^ - 







2x10"^ 




0x1 0""° 




-2x10'^ - 




10 12 14 16 18 20 
R (kpc) 

Figure 9. Hydrogen flow as function of R predicted by our 
best model with condensation of the corona. Inflows (outflows) 
are represented by positive (negative) values. Top panel: net flow 
(inflow — outflow) produced by the Galactic fountain. Bottom 
panel: accretion of pristine gas from the corona onto the disc. 
The global accretion rate is l.GMQyr"-*^. 



4.3 Accretion and circulation of gas 

Infall of metal-poor gas to the star-forming disc is an essen- 
tial ingredient of current models of the Galaxy's chemical 
evolu tion (e.g. iChiappini et al.. i2001; S chonri ch & Binne^ 
I2OO9I ). The predictions of such models depend to a signif- 
icant extent on the radial profile of the infall, but hitherto 
there has been no credible way of determining this profile. 
The prediction of our model for this profile is shown in Fig. [9] 
(bottom panel). The curve shows the pristine gas that, con- 
densing from the corona onto the fountain cloud wakes, fol- 
lows the cloud orbits back to the disc. The shape of the 
accretion profile is due to the variation of the mass outflow 
and the orbital time with radius. At _R < 3 kpc the specific 
accretion rate essentially vanishes because the orbits of the 
fountain clouds are confined within few hundred of parsecs 
from the disc (see Fig.[Sl). It then rises to a peak at i? ~ 9 kpc 
as in this region both the orbital time and the star formation 
are sufficiently high. At i? > 9 kpc the accretion rate falls 
again due to the low level of star formation at these radii, 
and drops to about zero beyond R = 13 kpc. The integral 
of this curve over the disc surface gives a global accretion 
rate of M = 1.6M0yr~^. Note that the peak in the accre- 
tion rate lies well beyond the peak in the Milky Way's SFR, 
which occurs around _R = 4kpc. This implies the need for a 
redistribution of gas in the disc, the study of which is be- 
yond the scope of the prese nt paper (but see for in stance 
iSchonrich fc Binnev 2009; Sp itoni fc Matteuccillioill '). 

The top panel of Fig. |9]shows the net flow (H I-|-H II) as 
function of R, obtained as the difference between the rate 
at which gas arrives at a given radius and the rate at which 
supernovae eject gas from that radius. The net flow profile 
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difTers from the accretion profile because fountain clouds do 
not fall back onto the disc at the same radius they are ejected 
from. As R increases, the orbits of fountain clouds start to 
get an inward component because of the interaction with 
the corona, thus clouds land at smaller and smaller radii. In 
particular, most of the clouds ejected at 11 < 7? < 15 kpc fall 
back onto the disc at 6 < i? < 11 kpc. As a consequence, a 
'draining region' forms at 11 <i?< 15 kpc, while inside this 
annulus the amount of gas falling onto the disc significantly 
exceeds the outflowing counterpart: about 70% of this inflow 
comes from the condensation of the corona, the remaining 
part is due to the described gas circulation. On the whole, 
the fountain circulation contributes to 76% of the halo mass, 
while the accreted gas is 24%. This latter is similar to that 
derived by FB08 for NGC 891 and NGC 2403 (~ 10 - 20%). 

4.4 HVCs and IVCs 

At large line-of-sight velocities the HI emission in the halo 
region of our Galaxy is dominated by the HVCs, the emission 
from which lies in rather isolated outlyi ng regions of (Z, b, v) 
space delimited by |wdev| > 90kms~^ (Ivan Woerden et al.l 
|200j). The upper panels in Fig. [10] compare the emis- 
sion predicted by our best model (left panel) near ulsr = 
— 144kms~^ with what is actually observed. In the observa- 
tional data, three main islands of emission are apparent. At 
b > 30° we see Complexes A and C, which have no real coun- 
terparts in the model. Similarly emission at I > 180° associ- 
ated with the Anticentre High- Velocity Cloud (ACHV) has 
no counterpart in the model. It is well known that the metal- 
licity of most HVCs, in particular Complexes A and C is 
much lower than that of gas in the disc, which stro ngly sug- 
gests that they come from outside the Galaxy (see I Wakkeil 
I2OOII ). Moreover, ass uming a distance of 10 kpc for both 
Com plexes A and C (jvan Woerden et~alll2004l : pThom et al.l 
l2008t ) gives for these objects heights from the midplane of 
6 — 8 kpc. To produce emission so far from the plane, we 
would need a kick velocity hv > 150kms~^. A model with 
hv so large completely fails to reproduce the data glob- 
ally. In summary, HVCs such as Complexes A and C are 
almost certainly extragalactic in origin and our model is 
quite correct not to reproduce them. The rate at which the 
Galaxy accretes ga s from infalling HVCs is ^ 0.2 Mqjv'^ 
ijSancisi et alj|200a ). an order of magnitude lower than our 
estimate of the rate of accretion via the fountain -I- corona 
condensation mechanism. 

The main body of emission in the upper right panel 
of Fig. \TU\ forms a large island that straddles the plane. Its 
summit lies at {l,b) = (120°, 15°) and is labelled "Outer 
Arm" . This feature is thought to arise from the warp of the 
HI disc. Our model (left panel) does not include a warp 
in the disc, so apart from two horns of very faint emission 
either side of the plane, its emission is confined in a thin 
region around the plane. 

The lower panels of Fig. [10] show predicted (left) and 
observed (right) emission at much lower heliocentric veloc- 
ities (wlsr — — 66kms~^) than the upper panel. At such 
Intermediate Velocities (IVCs are classically defined to be 
at 35kms~^ < I^devI < 90kms^^) the model is seen to 
reproduce the global structure of the data very well. More- 
over, some of the IVCs that contribute to this emission have 
been shown to have a distance from the Sun of < 3 kpc and 



disc- like metallicities l|van Woerden et al.ll200'3 ). consistent 
with their being fountain clouds. In the lower-right panel of 
Fig. \W\ several IVCs are visible (IV-Arch, IV-Spur, LLIV- 
Arch, AC-Shell, Complex K), mainly at positive latitudes. 
The emission of our model is smoother than the data, but 
apart from that it is consistent with the existence of all the 
IVCs. 

The classical IVCs shown in the lower right panel of 
Fig. \TU\ must be just the most conspicuous members of a 
large population of IVCs that collectively comprise the H I 
halo (see also MFll). By assuming azimuthal symmetry, 
our model suppresses much of the noise inherent in observ- 
ing a population of discrete clouds. In principle one could 
hope with a more sophisticated model to reproduce the sta- 
tistical properties of this noise, but there is no particular 
merit in reproducing individual IVCs, which will just be the 
chance, and ephemeral, products of individual bursts of star 
formation. 

If IVCs contain a mixture of metal-rich gas ejected 
from the disc and gas condensed from the rather metal-poor 
corona, the metallicities of these clouds should be intermedi- 
ate between those of the disc and corona. Since the fraction 
of accreted material is ~ 24% of the whole halo mass, the 
IVCs should be somewhat more metal-poor than the disc. 



5 DISCUSSION 

5.1 Analytic approximation vs hydrodynamical 
simulations 

In Section 12.31 we pointed out that our assumption of ex- 
ponential growth of the mass o f the cold gas is qual ita- 
tively similar to that observed by iMarinacci et all (|2010D m 
their hydrodynamical simulations of cold clouds travelling 
through a hot medium. Here we show that this agreement 
is quantitatively sound. 

The points in the upper panel of Fig. [TT]show the mass 
of coronal gas tha t has T < 3 x 10'* K in the simulations 
of iMarinacci et al] (|201ll ): the units are fractions of the ini- 
tial mass of the cloud. The curve shows the mass accreted 
by a particle that after 38 Myr starts to accrete accord- 
ing to equation ((SJ with the growth rate set to the value, 
a = 6.3 Gyr~^ , determined by our fits to the LAB datacube 
(Table [TJ. T he fit in Fig, j lllbetween the curve and the data 
points from iMarinacci et al.i (i2O10D is excellent. We inter- 
pret the delay by 38 Myr before accretion starts as the time 
required for cold gas from the cloud to mix with the coro- 
nal gas plus the time required for the coronal gas to cool to 
T < 3 X 10^ K. 

The points in the lower panel of Fig. [11] show the 
velocity of the cold-gas centroid in the simulations of 
iMarinacci et al.l (|201lf ). The dashed line shows the evolu- 
tion of the velocity if the cloud were a rigid body that ex- 
perienced hydrodynamical drag as describ ed in FB06 with 
^drag = 800 Myr. This line fits the data from iMarinacci et al.i 
(I2OIII ) well until t ~ 20 Myr. At later times the fuU thick 
curve shows the prediction for the combined effects of drag 
and condensation, obtained from equation ([7]l using again 
a — 6.3Gyr~^. The required condensation rate is very sim- 
ilar than that predicted by the simulations. 

It is remarkable that detailed hydrodynamical simula- 
tions yield a value of the accretion rate a that agrees so 
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Figure 10. Channel maps at j;lsr = ~144km (top) and tiLgj^ = — 66 km (bottom) for our best model with coronal condensation 
(left panels) and the LAB data (right panels). Both model and data are smoothed to 4° resolution. Contour levels in brightness 
temperature range from 0.04 K to 81.92 K scaling by a factor 2. Some of the high and intermediate-velocity complexes are labelled (see 
text). 



well with the value we obtained by fitting the H I datacube 
because the only connection between these two determina- 
tions is the underlying physics of turbulence. It is true that 
to obtain the agreement we did have to choose a value for 
the delay between the start of a hydrodynamical simulation 
and the onset of effective cooling. Some delay is inevitable 
because the hydrodynamical simulations start with unreal- 
istically spherical clouds and no wake. More realistic initial 
conditions for the dynamics of gas expelled from the disc by 
a superbubble would yield a time shorter than 30 — 40 Myr 
for steady accretion to become established. 

The simulations of iMarinacci et ahl l|201ll ) last only 60 
Myr; new simulations that also include gravity are now being 
carried out (Marinacci et al., in preparation). These show 
that the cold gas (cloud + wake) falls back to the disc in 
rather coherent structures with sizes of 1 — 1.5 kpc and after 
times also larger than 60 Myr and comparable to the ballistic 
travel times. 

While our choice of temperature threshold at T = 
3 X 10* K for gas to become visible as H I is somewhat arbi- 
trary, experimenting with different temperature thresholds 
we found that the trends of mass and velocity are perfectly 
comparable with the one shown here (F. Marinacci, private 
communication) . 

In conclusion, the analytical treatment of an inelastic 
collision between fountain particles and ambient gas pro- 



posed by FB08 and adopted her e is quantitatively supp orted 
by the numerical simulations of lMarinacci et al.l(|201ll ). The 
system of cloud -I- turbulent-wake accretes mass and looses 
specific momentum in the way predicted by the analytical 
formula. Remarkably the rate at which this process occurs 
in the simulations is the same as the rate we estimated by 
fitting models to the LAB survey data. 

5.2 Comparison to other galaxies 

FB06 and FB08 applied the model discussed here to two ex- 
ternal galaxies, the Milky Way-like NGC891, and the M33- 
like NGC 2403. For the former the kick velocities necessary 
to reproduce the data were in the range 70 — 90kms~^ de- 
pending on the shape of the potential and the extent and 
type of gas accretion considered. For NGC 2403 the data 
could be reproduced with values of the initial kick velocity 
from 75kms~^ down to about 50kms~^. Thus the value of 
70kms~^ that minimises the residuals in the Milky Way is 
fully consistent with those needed to produced the H I halos 
of other galaxies. As discussed in FB06, these velocities are 
also consistent with h ydrodynamical simula tions of super- 
bubble expansion (e.g. iMac Low et"allll989l ). 

The accretion parameters a estimated by FB08 for 
NGC 891 and NGC 2403 are a factor 3-7 lower than that 
found here for the Milky Way. This would imply shorter 
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Figure 11. Top panel: The square points show the mass of 
coronal gas that condensed into the cloud wake as a function of 
time extracted from the simulations of Marinacci et al. (2011). 
Bottom panel: The square points show the velocity centroid of 
cold gas in the simulations of Marinacci et al. (2011), the dashed 
curve shows the prediction in case of drag (tdrag = 800 Myr) 
without condensation. In each panel the solid curves show the 
predictions of the inelastic coUision+drag recipe used in this work 
(equation [7} ■ 



timescales for the condensation of the ambient gas in the 
Milky Way. However, the value of a is linked to the rota- 
tional velocity assumed for the corona where it interacts with 
the fountain. FB08 assumed that the fountain gas interacted 
with a static ambient medium, which maximised the effect 
of the process. Had we considered such a configuration for 
our corona we would have found a condensation parameter 
a = 2.5Gyr~^. In this case tdrag ^ 300, which implies that 
the drag dominates over the condensation. However, as ex- 
plained in Section [2.31 a static (non- rotating) corona is not 
a realistic possibility. 



6 CONCLUSIONS 

Just as in the 1980s photometry of external galaxies proved 
the key to making sense of star counts within our Galaxy 
l|Bahcall fc Soneiralll980D . so it is natural to turn to studies 
of the HI distributions of external galaxies for help in in- 
terpreting the Galaxy's HI datacube. In this paper we have 



investigated the extent to which the latter can be under- 
stood using a model of the extraplanar HI that emerged 
from observations of external galaxies such as NGC 891 and 
NGC 2403. In this model clouds of relatively cool gas are 
ejected from the plane by supernova-driven superbubbles 
and subsequently orbit over the disc on near ballistic trajec- 
tories whilst weakly interacting with the virial-temperature 
coronal gas through which they move. The interaction with 
the coronal gas is through a combination of ram-pressure 
drag and accretion of gas that cools in the cloud's wake. As 
a result of this accretion, more gas returns to the plane than 
left it. 

In light of recent developments, we modified the model 
slightly before applying it to the Galaxy. The principal mod- 
ifications were (a) to allow for rotation of the corona, and (b) 
to change the rule for determining the SFR from one based 
on total gas density to one based on the density of molecu- 
lar gas alone. Neither modification arises from a desire to fit 
data for the Galaxy better; (a) refiects an improved under- 
standing of the hydrodynamics of cloud-corona interaction, 
and (b) arises from better data for external galaxies. A task 
for the future is to reanalyse the data for external galaxies 
using our present model. 

Our model has three adjustable parameters, the char- 
acteristic velocity of cloud ejection h^, the specific accretion 
rate of clouds a, and a dimensionless parameter /ion, which 
determines where along its trajectory the cloud's gas, which 
is initially photoionised, becomes visible as HI. To fit the 
data we require ~ 70kms~^, regardless of the assump- 
tions we make about the values of a and /ion- 

At high Galactic latitudes, the observations show a clear 
bias towards negative velocities, and from this fact it follows 
that /ion has to be larger than zero. If we set a = 0, thus 
excluding accretion from the corona, we obtain a reasonable 
fit to the data with /ion = 1, so clouds become neutral only 
when they start to move back towards the plane. 

Permitting a to be non-zero yields fits to the data that 
are improved in small but significant ways. In particular, 
with a > 0, more emission is predicted at negative veloci- 
ties at low latitudes and / = 180°. Also less emission is pre- 
dicted at negative velocities and \l\ ~ ±10° and |6| ~ 30°. 
Both these changes arise from an increase in the extent of 
global infiow of the H I, and they improve the fit to the data. 
Moreover, with a the optimum value of /on is reduced 
from unity to 0.3, so clouds become neutral about 30 per- 
cent of the way to their highest point above the plane. Our 
best-fitting model has a = 6.3Gyr~^, which implies that 
1.6MQyr~^ of hydrogen is accreted from the corona. In- 
cluding the He content, this value rises to 2.3 Mq yr~^. This 
accretion rate is in excellent agreement with estimates of the 
accretion rate required to sustain the Galaxy's current rate 
of star formation without depleting its rather meagre stock 
of interstellar gas. 

Unfortunately, two problems prevent us from tightly 
constraining the value of a. The first is that it is inferred 
from some quite subtle features in the HI datacube. The 
second is that the optimum value of a depends on the 
value adopted for uiag, the equilibrium difference in the 
rotation velocities of the HI halo and the corona. When 
wiag is raised to lOOkms"^, the optimum value of a de- 
creases to 4.0 Gyr~^, which corresponds to an accretion rate 
of 1.0 Mq yr"^ (1.4 Mq yr"^ including the He content). 



14 A. Marasco, F. Fraternali, & J. J. Binney 



Fig. [To] illustrates how successfully the model simulates 
HI emission at Intermediate Velocities - the simulation is 
probably as perfect as it can be without reproducing individ- 
ual superbubbles. By contrast, the model does not reproduce 
emission at High Velocities, presumably because HVCs are 
extragalactic in origin. 

We find remarkable agreement between the optimum 
values of the model's parameters and the values found in 
earlier work. In particular, our value of ftv lies within the 
ranges of values found for NGC 891 and NGC 2403, and al- 
though our value for a is higher than that found for NGC 891 
this difference can be traced to different assumptions re- 
garding the rotation of the corona here and in the work of 
FB08 on NGC 891. Our value of a does beautifully repro- 
duce the results of a simulation by iMarinacci et al.1 (|201ll ) 
of the hydrodynamics of cloud-corona interaction. Our value 
dv4,/dz ~ — 14.3kms~^ kpc~^ for the vertical gradient of 
the halo's mean-streaming velocity agrees perfectly with the 
value (~ -15 kms"^ kpc"^) determined for NGC 891. These 
quantitative agreements between studies that use either rad- 
ically different methodologies or data inspires confidence 
that our underlying physical picture is correct. 

If we accept the fundamental soundness of the model, 
we can from it infer the three-dimensional structure of the 
HI halo, which is otherwise unknown. We find that the halo 
contains 3.0 x 10** Mq of gas, of which 2.7 x 10* Mq is 
neutral, in agreement with the value of 3.2^|!| g x 10* Mq es- 
timated by MFll. Its vertical scaleheight increases roughly 
linearly with radius from ~ 250 pc at J? = 3 kpc to 2 kpc at 
R = 14 kpc. At larger radii the scaleheight falls along with 
the local SFR. This HI halo is neither in hydrostatic equi- 
librium nor strictly in circular rotation: an analysis of the 
data that is nevertheless based on the assumption of hydro- 
static equilibrium and circular motion will yield a misleading 
structure for the HI disc, especially the flare. This structure 
will in turn lead to false conclusions about the distribution 
of matter required to gravitationally confine the H I disc (see 
Sect. 14.2]) . The rate per unit area at which the halo deposits 
coronal gas on the disc increases roughly linearly from near 
zero at i? = 3 kpc to a peak at R = 9 kpc and then falls 
roughly linearly to near zero at R = 13 kpc. Models of the 
chemical evolution of the disc need the accretion profile as 
an input, and it will be interesting to discover how the pre- 
dictions of such models change when the present accretion 
profile is adopted. 
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